COMPUTER METHOD AND SYSTEM FOR DETERMINING STRUCTURAL INFORMATION ON AN UNDERGROUND FORMATION AND CO

专利摘要:
Methods and systems for reconstructing low frequency particle velocity wave fields and eliminating ghosting in seismic cable data The present invention relates to methods and computational systems for eliminating ghosting in marine seismic cable data. in particular, an exploration seismology vessel tows numerous seismic cables that form a data acquisition surface located below a free surface. the methods computationally remove phantom or substantially remove phantom receiver signals from the seismic data recorded by the seismic cable receivers. ghosting elimination methods include low frequency compensation to recover vertical velocity wave field information that is normally lost due to a low signal-to-noise ratio over a low frequency range regardless of free surface conditions or shape of the surface. data acquisition surface. .
公开号:BR102013017551B1
申请号:R102013017551-0
申请日:2013-07-09
公开日:2021-08-03
发明作者:Stian Hegna；Walter Söllner；Okwudili Orji
申请人:Pgs Geophysical As；
IPC主号:

专利说明:

BACKGROUND
[0001] The present invention relates to, in recent decades, the petroleum industry having invested heavily in the development of marine seismic survey techniques that provide knowledge of underground formations under a body of water to seek and extract valuable mineral resources such as like oil. High resolution seismic images of an underground formation are essential for quantitative seismic interpretation and improved reservoir monitoring. For a typical marine seismic survey, an exploration seismology vessel tows one or more seismic sources and one or more cables below the water surface and over an underground formation to be analyzed for mineral deposits. The ship contains seismic acquisition equipment such as navigation control, seismic source control, seismic receiver control and recording equipment. Seismic source control causes one or more seismic sources, which are normally air guns, to produce acoustic impulses at selected times. Each impulse is a sound wave that travels through water and underground formation. At each interface between different rock types, a portion of the sound wave is refracted, a portion of the sound wave is transmitted and another portion is reflected back towards the water body to propagate towards the surface. The towed cables behind the ship are elongated cable-like structures. Each cable includes numerous seismic receivers or sensors that detect pressure and/or changes in particle movement in the water created by sound waves reflected back into the water from underground formation.
[0002] Sound waves that propagate upwards from underground formation are referred to as "upward going" wave fields that are detected by receivers and converted into seismic signals that are recorded by recording equipment and processed to produce images seismic characteristics that characterize the geological structure and properties of the underground formation being analyzed. However, seismic signals can also include "source ghosting" produced by sound waves that are first reflected off the sea surface before the waves travel on the subsurface to produce scattered wave fields detected by receivers. Source ghosts are time-delayed relative to sound waves that travel directly from the source to the underground formation. As a result, source ghosts can amplify some frequencies and attenuate other frequencies and typically manifest as spectral discontinuities in recorded seismic waveforms, making it difficult to obtain accurate high-resolution seismic images of underground formation. In addition to "source ghosts," the seismic signal can also include "receiver ghosts" produced by scattered sound waves that are first reflected off the sea surface before reaching the receivers. Receiver ghosting can also amplify some frequencies and attenuate other frequencies and are commonly manifested as receiver ghosting discontinuities. As a result, those working in the petroleum industry continue to look for systems and methods to remove the effects of ghost reflections or eliminate the ghosting of seismic signals. DESCRIPTION OF DRAWINGS
[0003] Figure 1 shows a domain volume of the earth's surface.
[0004] Figure 2 shows subsurface features of an underground formation in the lower portion of the domain volume shown in Figure 1.
[0005] Figures 3A to 3C show a method of exploration seismology whereby digitally encoded data is acquired instrumentally for subsequent exploration seismology processing and analysis to characterize the structures and distributions of features and materials that support the solid surface of the earth.
[0006] Figures 4A to 4C show processed waveforms generated from hydrophone and geophone outputs.
[0007] Figure 5 shows an elevation view of an exploration seismology vessel that tows a source and a cable located below a free surface.
[0008] Figure 6 shows a flowchart associated with a method for eliminating the ghosting effect of seismic signal data source.
[0009] Figures 7A and 7B introduce coordinates and terminology associated with a computational process for eliminating the ghosting of the font.
[00010] Figure 8 shows an elevation view of a seismology exploration vessel and a cable and includes an enlarged view of a receiver located along the cable.
[00011] Figures 9A to 9C each show an exemplary plot that represents an aspect of wave field decomposition.
[00012] Figures 10A to 10C each show an exemplary plot that represents an aspect of computing an imaged free surface profile.
[00013] Figure 11 shows a control flowchart of a computational method to calculate a free surface profile.
[00014] Figure 12A shows a plot of a hypothetical free surface and an imaged free surface profile that approximates the free surface.
[00015] Figure 12B shows a plot of a hypothetical free surface and an extended imaged free surface profile that approximates the free surface.
[00016] Figure 13 shows parameters of an equation representing free surface wave field reflectivity.
[00017] Figures 14A to 14B show a two-dimensional view of an exemplary model space associated with a fluid volume.
[00018] Figure 15 shows a segment of a cable located below an imaged free surface of a fluid volume.
[00019] Figure 16 shows an example of a generalized computer system that performs an efficient method for eliminating the ghosting of a scattered wave field.
[00020] Figures 17A to 17E show plots of a cable and a free surface. DETAILED DESCRIPTION
[00021] Methods and computational systems for eliminating ghosting in marine seismic cable data receiver are described. In particular, an exploration seismology vessel tows numerous cables that form a data acquisition surface located below an air/fluid surface referred to as the "free surface." cables include receivers that measure pressure and particle motion wave fields that are digitally encoded and stored. The methods computationally eliminate ghosting or substantially remove receiver phantom signals from the seismic data recorded by the receivers regardless of free surface conditions or the shape of the data acquisition surface. In other words, the methods described below computationally remove receiver phantom signals from the seismic data without assuming restrictions on the free surface or assuming restrictions on the shape of the data acquisition surface, such as assuming a "frozen" flat free surface (ie, static ) and assume a "frozen" flat, horizontal data acquisition surface. Ghost elimination methods include low frequency compensation to recover vertical velocity wave field information that is normally lost due to low signal-to-noise ratio over a low frequency range.
[00022] The following discussion includes two subsections: in subsection (I), an overview of exploration seismology is provided; and in subsection (II), a discussion of computational processing methods for eliminating ghosting in seismic signal data receiver as an example of computational processing methods and systems to which this disclosure is directed. The reading of the first subsection may be omitted by those already familiar with marine exploration seismology. I. An overview of marine exploration seismology
[00023] Figure 1 shows a domain volume of the earth's surface. Domain volume 102 comprises a solid volume of sediment and stone 104 below the solid surface 106 of the earth which, in turn, lies beneath a fluid volume of water 108 within an ocean, an arm or bay or a large lake. of fresh water. The domain volume shown in Figure 1 represents an exemplary experimental domain for a class of observational and analytical exploration seismology techniques and systems referred to as "marine exploration seismology".
[00024] Figure 2 shows subsurface features of an underground formation in the lower portion of the domain volume shown in Figure 1. As shown in Figure 2, for the purpose of exploration seismology, the fluid volume 108 is relatively featureless, usually a homogeneous volume overlying the solid volume 104 of interest. However, while the fluid volume 108 can be explored, analyzed, and characterized with relative accuracy by using many different types of methods and probes, which include remote sensing submarines, sonar, and other such devices and methods, the volume of solid crust 104 underlying the fluid volume is comparatively far more difficult to probe and characterize. Unlike the overlying fluid volume 108, the solid volume 104 is significantly heterogeneous and anisotropic and includes many different types of features and materials of interest to exploration seismologists. For example, as shown in Figure 2, solid volume 104 may include a first layer of sediment 202, a first layer of fractured and raised stone 204, and an underlying second layer of stone 206 below the first layer of stone. In certain cases, the second layer of stone 206 may be porous and contain a significant concentration of liquid hydrocarbon 208 which is less dense than the material of the second layer of stone and therefore rises upwards within the second layer of stone 206 In the case shown in Figure 2, the first layer of stone 204 is non-porous and therefore forms a cover that prevents further upward migration of liquid hydrocarbon which therefore accumulates in a saturated hydrocarbon layer 208 below the first layer 204. An objective of exploration seismology is to identify the locations of porous hydrocarbon-saturated extract within volumes of the earth's crust underlying the earth's solid surface.
[00025] Figures 3A to 3C show a method of exploration seismology whereby digitally encoded data is acquired instrumentally for subsequent exploration seismology processing and analysis to characterize the structures and distributions of features and materials of an underground formation. Figure 3A shows an example of an exploration seismology vessel 302 equipped to perform a continuous series of exploration seismology experiments and data collections. In particular, the ship 302 tows one or more cables 304 to 305 through a plane of approximately constant depth located generally several meters below the free surface 306. cables 304 to 305 are long cables that contain power transmission lines and data to which receivers, also referred to as "sensors", are connected at regular intervals. In one type of scanning seismology, each receiver, such as the receiver represented by shaded disk 308 in Figure 3A, comprises a pair of seismic receivers that include a geophone that detects vertical displacement within the fluid medium over time by detecting a particle movement, velocities or accelerations and a hydrophone that detects changes in pressure over time. Cables 304 through 305 and vessel 302 include sophisticated electronic sensor equipment and data processing facilities that allow receiver readings to be correlated to absolute free surface positions and absolute three-dimensional positions relative to an arbitrary three-dimensional coordinate system. In Figure 3A, receivers along the cables are shown as at rest below free surface 306, with receiver positions correlated to positions on the overlying surface, such as a surface position 310 correlated to receiver position 308. also tows one or more acoustic wave sources 312 that produce pressure pulses at spatial and temporal intervals as ship 302 and towed cables 304 to 305 move along free surface 306. Sources 312 may also be towed by other ships or they can be otherwise arranged in the fluid volume 108.
[00026] Figure 3B shows an expanding spherical acoustic wavefront represented by semicircles of increasing radius centered on the acoustic source 312, such as a semicircle 316, following an acoustic pulse emitted by the acoustic source 312. The wavefronts are, in practice, shown in cross-section of the vertical plane in Figure 3B. As shown in Figure 3C, the downwardly expanding acoustic wave field shown in Figure 3B eventually reaches the solid surface 106, at which point the downwardly and outwardly expanding acoustic waves partially reflect from the solid surface and partially refract down into the solid volume, becoming elastic waves within the solid volume. In other words, in fluid volume, waves are compressive pressure waves or P waves, the propagation of which can be modeled by the acoustic wave equation whereas, in a solid volume, waves include both P waves and transverse waves or S waves, the propagation of which can be modeled by the elastic wave equation. Within the solid volume, at each interface between different types of materials or at discontinuities in density or in one or more of several other physical characteristics or parameters, downwardly propagating waves are partially reflected and partially refracted as a solid surface 106. As a result, each point on the solid surface and within the underlying solid volume 104 becomes a potential secondary point source from which elastic and acoustic waves, respectively, can emanate upward toward receptors in response to the emitted pressure impulse. by the acoustic source 312 and downwardly propagating elastic waves generated from the pressure impulse.
[00027] As shown in Figure 3C, secondary waves of significant amplitude are emitted generally from points on or near the solid surface 106 such as point 320 and from points at or very near a discontinuity in the solid volume 104, such as as points 322 and 324. Tertiary waves can be emitted from the free surface 306 back towards the solid surface 106 in response to secondary waves emitted from the solid surface and subsurface features.
[00028] Figure 3C also shows the fact that secondary waves are generally emitted at different times within a time range following the initial pressure impulse. A point on solid surface 106 such as point 320 receives a pressure perturbation that corresponds to the initial pressure pulse faster than a point within solid volume 104 such as points 322 and 324. Similarly, a point on solid surface directly underlying the acoustic source receives the pressure impulse faster than a point that lies farther away on the solid surface. Therefore, the times at which secondary and higher order waves are emitted from various points within the solid volume are related to the distance, in three-dimensional space, of the points from the acoustic source.
[00029] Elastic and acoustic waves, however, travel at different speeds within different materials as well as within the same material under different pressures. Therefore, the travel times of the initial pressure pulse and secondary waves emitted in response to the initial pressure pulse are complex functions of the distance from the acoustic source as are the materials and physical characteristics of the materials by which the acoustic wave corresponding to the pressure pulse initial travels. Furthermore, as shown in Figure 3C for the secondary wave emitted from point 322, the shapes of the expanding wavefronts can change as the wavefronts cross interfaces and as the speed of sound varies in the media traversed by the wave. The superposition of waves emitted from within domain volume 102 in response to the initial pressure pulse is a generally very complicated wavefield that includes information about the shapes, sizes and material characteristics of domain volume 102 that includes information about the shapes, sizes and locations of the various reflective features within the subterranean formation of interest to exploration seismologists.
[00030] The complicated wave field that ensues in response to the initial pressure impulse is sampled over time by receivers positioned along the cables towed by a seismology exploration vessel. Figures 4A to 4B show processed waveforms generated by a hydrophone and a geophone, respectively. As shown in Figure 4A, the waveform recorded by the hydrophone represents the pressure at moments following the initial pressure pulse with the amplitude of the waveform at a point in time related to the pressure on the hydrophone at the point in time. Similarly, as shown in Figure 4B, the geophone provides an indication of particle motion or fluid velocity or acceleration in a vertical direction with respect to time. The pressure and particle motion signals represented by the waveforms in Figures 4A and 4B can be used to generate wave fields from pressure and vertical velocity wave fields which, in turn, can be used for seismic processing, such as as attenuation of multiples in marine seismic data. However, because recorded particle motion signals are often contaminated with low frequency noise due to vibrations from towed cables, the signal-to-noise ratio over a low frequency range for the combined signals is poor. Figure 4C shows a plot of relative amplitude against a range of acoustic frequencies of the exemplary pressure and particle motion signals depicted in Figures 4A and 4B. Solid curve 401 and dotted curve 402 represent the recorded pressure and motion signals, respectively. Figure 4C depicts a portion 403 of the particle motion signals with relative amplitudes greater than zero, which correspond to a low signal-to-noise ratio over a low frequency range of 0 Hz at an acoustic threshold frequency
which in the example in Figure 4C is about 20Hz. As a result, the particle movement signal below the threshold frequency
it cannot normally be used to perform seismic processing. On the other hand, the pressure wave field below the threshold frequency
can be used to recreate certain information lost by the particle motion sensor when the pressure signal spectrum has a satisfactory signal-to-noise ratio and when the depth of the pressure and particle motion sensors is known. As a result, the cable depth is selected so that the frequency of the first spectral discontinuity in the pressure sensor signal caused by surface reflections is greater than the threshold frequency.

[00031] The free surface acts as an almost perfect acoustic reflector, which causes ghosting effects on the recorded seismic data. At the location of the source of each pulse, a time-delayed reflection from the free surface, called a "source ghost," chases the seismic wave field that travels directly from the source into the earth's subsurface. As a result, both low and high frequency information are penalized and the earth's subsurface cannot be accurately imaged at all seismic frequencies. In addition, at each receiver location along a cable, a time-delayed reflection of the free surface, called a "receiver ghost," interferes in a continuous and unwanted manner with the seismic wave field spread directly from the subsurface of the ground to the receiver. Figure 5 shows an elevation view of an exploration seismology vessel 502 towing a source 504 and a cable 506 located below a free surface 508. Figure 5 includes four ray paths representing four ways in which wave fields associated with the same pulse may travel from source 504 to a receiver 510 located along cable 506. The dashed line 514 represents a primary ray or direct path in which a pulse output from source 504 travels directly on subsurface 512 and fields of Associated scattered waveforms travel directly from the subsurface 512 to the receiver 510. The solid line 516 represents a ray path that includes a source ghost 517 produced by sound waves that are first reflected from the free surface 508 before entering the subsurface 512 to produce scattered wave fields that travel directly from the subsurface 512 to the receiver 510. The long dashed line 51 8 depicts a ray path including a receiver ghost 519 produced by sound waves traveling directly from source 504 to subsurface 512 to produce scattered wave fields that are reflected from free surface 508 before being detected by receiver 510. The dashed dotted line 520 represents a ray path that includes both a source phantom 521 and a receiver phantom 522.
[00032] The scattered wave fields associated with the primary ray paths are ideally the desired seismic wave fields to be detected by the receivers. In practice, however, source phantom and receiver phantom are also detected and because phantoms are delayed in time, phantoms can constructively and destructively interfere with the recorded waveforms of the seismic wave field. As a result, ghosts can lead to inaccurate seismic images of underground formation located below the fluid volume. The methods and computational systems described below are aimed at eliminating ghosting in a seismic signal data receiver for an arbitrarily free surface that varies in time and without restrictions on the shape of the data acquisition surface defined by the cables. The methods are also not limited to using pressure wave and vertical velocity wave fields above the threshold frequency.
In other words, receiver ghosting elimination in seismic signal data can be achieved for pressure wave and vertical velocity wave fields above the full frequency range. II. Methods for Eliminating the Ghost Effect in Receptor as an Example of Computational Processing Methods and Systems
[00033] Figure 6 shows a control flowchart 600 associated with a computational method for eliminating the ghosting in a seismic signal data receiver with low frequency compensation. Flowchart 600 summarizes a general computational process for eliminating ghosting in a seismic data receiver. Each block of the 600 diagram is described below in a separate subsection. Estimate Cable Depth
[00034] In block 601 of Figure 6, the cable depth is estimated. Figures 7A and 7B introduce coordinates and terminology associated with a computational process for eliminating the ghosting in receiver. Figure 7A shows a top or xy plane view of an example exploration seismology vessel 702 that tows a source 704 and four separate cables 706 to 709 located below a free surface. Each cable being attached at one end to ship 702 and at the opposite end to a buoy such as a buoy 710 attached to cable 708. Cables 706 to 709 ideally form a horizontal planar acquisition surface located below the free surface. However, in practice, the acquisition surface can vary smoothly in the z direction due to ocean currents and active environmental conditions. In other words, towed cables can also curl as a result of dynamic fluid conditions. Figure 7B shows an elevation or xz plan view of cable 708 located below a free surface 712. Figure 7B represents a snapshot, at one instant in time, of the undulating free surface 712 and the corresponding smooth wave-like shape in cable 708. Figure 7A includes an xy plane 714 and Figure 7B includes an xz plane of the same Cartesian coordinate system used to specify coordinate locations within the fluid volume with respect to three spatially orthogonal geometric axes labeled x, y, and z . The x coordinate uniquely specifies the position of a point in a direction parallel to the length of the cables and the y coordinate uniquely specifies the position of a point in a direction perpendicular to the x axis and substantially parallel to the free surface 712 and the z coordinate uniquely specifies the position of a point perpendicular to the xy plane. As shown in Figure 7B, cable 708 is at a depth zr, which can be estimated at various locations along the cables from hydrostatic pressure measurements made by depth controllers (not shown) attached to the cables. Depth controllers are normally placed at intervals of about 300 meters along each cable. The estimated cable depths are then used to calculate a two-dimensional interpolated cable shape that approximates the wave-like shape of a real cable at one instant in time. Alternatively, the estimated cable depths can be used to calculate a three-dimensional interpolated surface approximation of the acquisition surface. A zr depth and elevation of the free surface profile are estimated relative to the geoid, which is represented in Figure 7B by the dotted line 718. The geoid is the hypothetical land surface that coincides anywhere with mean sea level and is used to define zero elevation (that is, z = 0). Shaded disks, such as shaded disk 720, represent receivers spaced at regular intervals - . The coordinates of receiver 720 are given by - , where a depth Zr can be an interpolated value. Measurement of Pressure and Velocity Wave Fields
[00035] Returning to Figure 6, in block 602, pressure and velocity wave fields are measured at the receivers. Figure 8 also shows an elevation view of ship 702 and cable 708 and includes an enlarged view 802 of a receiver 804 located along cable 708. Magnified view 802 describes that the receivers are actually dual sensors that include a sensor pressure 806 such as a hydrophone and a motion sensor 808 such as a geophone. Each pressure sensor measures a pressure wave field denoted by p(xr,yr,zr,t) and each of the motion sensors measures a velocity wave field denoted by
Where
and y represents the spatial coordinates x and y of the sensor and t represents time. The m index is a positive integer used to identify the receiver and is also referred to as a "channel" index. The acquisition surface, as shown in Figure 7A, has 4 cables 706 to 709 and each cable includes 14 receivers. Receivers can be indexed from 0 to 13. For example, a receiver 810 on cable 708 closest to ship 702 can be assigned a channel index m = 0 and a receiver 812 located along the same cable 708, farther away of ship 702, channel index m = 13 can be assigned with dual sensors indexed from 1 to 12 consecutively. Motion sensors can be mounted within a rocker arm to guide motion sensors to detect particle movement in a direction normal to the cable as represented by unity normal vector 814 in close-up view 802. Therefore, motion sensors measure a field of velocity wave directed normal to cable that varies smoothly 708, where subscript vector 7 represents a unit normal vector that points in the direction 814 normal to cable 608 in the xz plane. However, for measured particle motion below the threshold frequency .
the signal-to-noise ratio is typically low as described above with reference to Figure 4C. As a result, vertical speeds below the threshold frequency do not
can usually be accurately determined. Methods include the low frequency compensation described below to calculate vertical velocity wave fields on cable receivers below the threshold frequency
. wave field decomposition
[00036] Returning to Figure 6, in block 603, the wave field decomposition is performed on the pressure and velocity wave fields for frequencies greater than the threshold frequency
described above with reference to Figure 4C. In the following description the y spatial component is ignored to simplify the description. Note that, in practice, the y spatial component is included. In other words, in the discussion that follows, the three spatial coordinates of the measured pressure wave field
are reduced to two spatial coordinates in the pressure wave field representation.
and the three spatial coordinates of the measured velocity wave field
rj are reduced to two spatial coordinates in the velocity wave field representation.
. Reduction to two spatial coordinates gives clear insight while preserving the main features of the method.
[00037] Pressure and velocity wave fields can be decomposed into vertical velocity and pressure components that go up and down. In Figure 8, directional arrow 816 represents the direction of an upward-going wavefield detected by a receiver 818 and the dashed-line directional arrow 820 represents the direction of an upward-going wavefield reflected from the surface. free 712 and detected by receiver 818. In other words, the pressure wave field
it is composed of a pressure component going up and a pressure component going down and the velocity wave field.
it is also composed of a vertical velocity component going up and a vertical velocity component going down. The downward going pressure wave field and the downward going vertical velocity wave field are receiver phantom wave fields and the upward going pressure wave field and the going vertical velocity wave field upwards are wave fields with the receiver's ghosting effect.
[00038] Figures 9A to 9C each show an exemplary plot that represents an aspect of wave field decomposition. In Figures 9A-9C, vertical axis 902 is a z-coordinate geometric axis representing a depth or spatial dimension; and in Figures 9A-9B, horizontal axis 904 is an x-coordinate geometric axis representing a spatial dimension x . The spatial coordinate "0" z corresponds to a plane tangent to .■• = : to geoid 718 described above with reference to Figure 7B. Figure 9A shows an exemplary plot of a 906 interpolated cable based on the estimated depths obtained from depth controllers attached to the actual cable as described above with reference to Figure 7. The shape of the 906 interpolated cable substantially corresponds to the shape of the 708 cable shown in Figure 7B. Shaded disks, such as shaded disk 908, represent the locations of receivers along interpolated cable 906. The pressure wave field
the velocity wave field s
: associated with each receiver is transformed from the space-time domain to a pressure wave field
and a wave field of velocity r
in the space-frequency domain, where
is the angular frequency for the acoustic frequencies / detected by the dual sensors. For example, the pressure wave field and velocity wave field associated with a receiver can be transformed by using Fourier transforms:

[00039] where j is the imaginary unit
In practice, the transformation can be performed by using a Fast Fourier Transform ("FFT") for greater computational efficiency. Note that lowercase letters p and v are used to represent quantities in the spacetime domain while uppercase letters P and V are used to represent quantities in the space-frequency or Wavenumber-Frequency domain. Figure 9B shows an exemplary plot of the interpolated cable 906 with the pressure and velocity wave fields transformed to the space-frequency domain.
[00040] After the pressure and velocity wave fields associated with each receiver have been transformed from the space-time domain to the space-frequency domain, the pressure wave fields and the velocity wave fields are combined to produce an upward component of pressure in the geoid (ie, z = 0) in the wavenumber-frequency domain. Pressure sensors alone cannot distinguish between the opposite polarity of the scattered seismic wave field from the underground formation and the time-delayed seismic wave field reflected downward from the sea surface (ie receiver phantom), but the information obtained from the particle motion sensors can be used to distinguish the polarity of the seismic wave field. A mathematical expression relating the pressure and velocity wave fields to a pressure component that goes above the pressure wave field in the geoid is given by:

[00041] where kz is the vertical wave number in the z direction given by:

[00042] with c being the velocity of sound in the fluid;
[00043] kx is the horizontal wave number in the spatial direction x
[00044] m is the dual sensor or channel index;
[00045] M is the total number of dual sensors located along the cable;
[00046] p is the density of the fluid;
[00047] z r,m is the interpolated depth of the cable in the m-th dual sensor;
[00048] nx is the x component of the normal vector

[00049] nz is the z component of the normal vector
and
[00050]
is the velocity wave field for angular frequencies greater than a threshold angular frequency (ie, .-is-
)
[00051] Similarly, the component of pressure that goes down in the wave-frequency number domain at z = 0 is calculated in a similar way by:

[00052] Note that the pressure wave field going up high and the pressure wave field going down low are computed from the pressure wave field and the velocity wave rest field.

[00053] Figure 9C shows an exemplary plot of the wave-frequency number domain. The horizontal axis 910 represents the geometric axis of the kx coordinate of the wavenumber. Note that the geometric axis of angular frequency w is perpendicular to the plane kxz, which is not represented in Figure 9C. Shaded disk 912 represents a point in the wave-frequency number domain located along the geometric axis of the kx coordinate 910. The point in the wave-frequency number domain 912 has pressure plane wave field components that go to up and down
and
916.
[00054] After the pressure wave field components going up and going down in the geoid have been determined in the wave-frequency number domain, the pressure wave field components are shifted to an observation level in a zobs depth between the geoid and the cape. In Figure 9C, dashed line 918 represents the geometric axis kx at an observation level located between the geoid and the cable indicated on the z axis by zr 920. A zobs depth of observation level 918 rests between the geoid and a cable depth 920. The pressure component that goes up at a point
921 over observation level 918 zobs is calculated
dust

[00055] Similarly, the component of pressure that goes down at the zobs 918 observation level is calculated by:

[00056] The vertical velocity component that goes up at the point
921 is calculated from the upward pressure component
per:

[00057] The vertical velocity component that goes down at the point
921 is calculated from the downward pressure component
per:

[00058] In other embodiments, the calculation of the free surface profile can also be achieved by using the vertical velocity components going up and down.
[00059] In short, wave field decomposition is the process of transforming the pressure wave field
and the velocity wave field
measured on the corresponding dual sensors of each receiver in the space-time domain, as shown in Figure 9A, in a pressure component that plowerj, 7obs goes down
922, an upward pressure component;
923, a vertical velocity component that goes down -
924 and a vertical velocity component going up -
925 in the wave-frequency number domain as shown in Figure 9C. The high, low, high, and low are computed from the pressure wave field P and the velocity wave field
above the threshold frequency
described above with reference to Figure 4C.
[00060] Note that a three-dimensional spatial version of the pressure components going up and down and the vertical velocity components going up and down can be obtained by substituting the vertical wavenumber

[00061] where ky is the horizontal wave number in the y spatial direction. Compute a Free Surface Profile
[00062] Returning to Figure 6, after wave field decomposition has been achieved in block 603, an imaged free surface profile is calculated in block 604. Figures 10A to 10C show three different plots in a profile computation process free surface imaged from a real free surface. An imaged free surface profile is first calculated by computing pressure components that go up and down or vertical velocity components that go up and down for each point in a series of depth levels that extend upwards. from the observation level beyond the geoid or geometric axis 910. Figure 10A shows an example of a series of points represented by a column of open circles 1002 that extend from point 921 on the zobs observation level 918 beyond of geometric axis 910. Each open circle in the 1002 series represents a point with coordinates , where the kx coordinate is the same for each point and the z coordinates grow incrementally in a process known as extrapolation. The upward pressure wave field is calculated at each point in the 1002 series by:

[00063] where z is a coordinate value in series of points 1002. The pressure wave field going down is calculated at each point in series 1002 by:

[00064] In another modality, the vertical velocity component that goes up at each point in the 1002 series is calculated by:

[00065] The vertical velocity component that goes down at each point in the 1002 series is calculated by:

[00066] After the extrapolated up and down pressure components or the extrapolated up and down vertical velocity components have been calculated in the wave-frequency number domain, an inverse Fourier transform is used to transform the extrapolated up and down pressure components extrapolated and/or extrapolated up and down vertical velocity components extrapolated in the space-frequency domain. Figure 10B shows a series of points 1010 in the space-frequency domain described above with reference to Figures 9A to 9B. An inverse Fourier transform is used to transform the pressure and velocity components associated with points in the 1002 series in Figure 10A to obtain corresponding pressure and velocity components of points 1010 in Figure 10B. The transformed pressure wave fields for the space-frequency domain are given by:

[00067] The transformed vertical velocity wave fields to the space-frequency domain are given by:


[00068] For example, an inverse Fourier transform can be used to transform the vertical velocity and pressure components that go up and down 1005 to 1008 associated with point (kx,z) 1004 shown in Figure 10A into a component of pressure going down plow(xr,z,w) 1013, a component of pressure going up . pcup(xr,z,w) 1014 a vertical velocity component that goes down vzdown(xr,z,w) 1015 and a vertical velocity component that goes up vzup(xr,z,w) 1016 at a corresponding point (xr,z) 1018 in the 1010 series shown in Figure 10B. In practice, the transformations represented by equations 11a through 11b and 12a through 12b can be performed using a discrete Inverse Fast Fourier Transform ("IFFT") for greater computational efficiency.
[00069] An imaging condition is used to calculate an image value I(xr,z) at each point in the series of points in the space-frequency domain. For example, in Figure 10C, an imaging condition is applied to each point in the series of points 1010 to obtain an image value such as an image value I(xr,z) associated with the point (xr,z)1018. The imaging condition can be a cross-correlation of vertical velocity or pressure components going up or down extrapolated in the space-frequency domain. In one embodiment, the imaging condition that represents a free surface image value for a selected receiver x position and z extrapolation depth is calculated by applying the following cross-correlation equation:

[00070] where the overlying bar designates a complex conjugation. In one modality,
represents
and
represents
. In another modality,
represents
represents
. In other modalities, the imaging condition can be a normalized cross-correlation given by:

[00071] The z coordinate value or extrapolation depth associated with the maximum image value. Imax(xr,z) for a given channel position x corresponds to a free surface elevation z at the position of the receiver x. In the example in Figure 10C, the points on the 1010 series associated with the 908 receiver are represented as dashed-line open circles except for point 1022 which corresponds to a point (xr,z) that yields maximum cross correlation.
. As a result, the point (xr,z0) is an image point of a free surface profile. Modalities are not limited to the imaging conditions described above. Other types of imaging conditions can also be used. Other image points represented by open circles such as 1024 open circles are calculated in the same way by applying the same imaging condition to each point in a series of points as described above. The interpolation over the collection of image points forms an approximate static time image of the free surface profile above cable 908 that corresponds to a selected time window of the input data. For example, in Figure 10C, a curve 1026 is a free surface profile imaged from the hypothetical free surface profile represented by the dotted curve 1028 in a selected time window. A free surface profile imaged above a cable can be interpolated by using spline interpolation, Lagrange interpolation, polynomial interpolation, or another suitable interpolation technique. In other embodiments, image points associated with two or more cables can be used to calculate an approximate three-dimensional free surface above the cables by using multidimensional interpolation techniques such as Barnes interpolation, Bezier interpolation, bicubic interpolation, and bilinear interpolation. The points along the free surface profile image 1026 are represented by [x,f(x)] where x represents a horizontal coordinate along the geometric axis x 904 and f(x) is the value of the interpolated function representing height of the free surface. A time-varying free surface estimation is obtained by continuously moving the input data time window and by combining the computed "frozen" free surface profiles for all time windows from a recording start time to a recording end time.
[00072] Figure 11 shows a control flowchart of a computational method to calculate the free surface profile of block 604. In block 1101, vertical pressures and velocities that go up and down are calculated for points in the number domain of wave-frequency as described above with reference to Figure 10A. In block 1102, the vertical pressures and velocities associated with the points are transformed from the wavenumber-frequency domain to the space-frequency domain. In block 1103, an imaging condition is applied as described above with reference to Figure 10C. In block 1104, the point associated with the maximum imaging condition for a given receiver channel is calculated, as described above with reference to Figure 10C. At block 1105, the image points associated with a cable are used to calculate a "frozen" image profile of the free surface. Free Surface Profile Extension
[00073] Returning to Figure 6, at block 605, the free surface profile is extended to cover the source location. Figure 12A shows a plot of a hypothetical free surface represented by a dashed curve 1028 and the imaged free surface profile 1026 that approximates the free surface 1028 above cable 906. Circle 1202 represents the coordinate location of a source towed by a exploration seismology vessel (not shown). As shown in Figure 12A and Figures 5, 7B and 8 described above, source 1202 is not located over the cables, but rather is located between the ship and the cables. Because the imaged free surface profile 1026 is computed from subsurface reflections, the imaged free surface profile 1026 is limited to cable scattering as indicated by dashed lines 1204 and 1206. As a result, the imaged free surface profile 1026 does not covers the 1202 source and can only be used to ghost out of ghosting that is reflected from the free surface directly above the source, and cannot be used to ghost out of ghosting that is reflected from the free surface above the 1202 source. For example, directional arrows 1208 and 1210 represent two outputs of source sound impulse paths 1208 and 1210. Both paths result in straight downward going wave fields with ghosting and source ghost reflections as described above with reference to Figure 5. The computation of the straight-down wave field that includes its phantom is used for field reconstruction of low frequency vertical velocity waveform and elimination of the ghosting of seismic data with sources above cable level. The imaged free surface 1026 can be used to eliminate the ghosting of seismic data with ghosts following path 1210, but because the imaged free surface 1026 does not extend to cover the region above the source 1202, the free surface image 1026 does not can be used to eliminate the ghosting of seismic data that travel along path 1208.
[00074] To compensate for ghosting associated with free surface reflections above the source, the imaged free surface profile is extended by calculating a free surface profile extent of the imaged free surface profile that covers the region above the coordinate location of the source. Figure 12B shows a plot of the hypothetical 1028 free surface and an extended imaged free surface profile. The extended imaged free surface profile is the imaged free surface 1026 extended to include a free surface profile extension 1212. The free surface profile extension 1212 can be calculated based on parameters associated with a realistic physical model of a state of the body. fully developed sea determined at the time the seismic data is measured. For example, a realistic physical model might be a Pierson-Moskowitz model that provides parameters such as wavelength, wave heights, and the speed of free surface waves. The Pierson-Moskowitz free surface model can be used to calculate the extent of free surface profile 1212 that covers source 1202.
[00075] The Pierson-Moskowitz model assumes that when the wind blows steadily for a long period of time over a large area, the waves eventually reach a state of equilibrium with the wind. This condition is referred to as a "fully developed sea". The Pierson-Moskowitz spatial roughness spectrum for a fully developed sea in one direction is given by:

[00076] where K is the space wave number of the free surface wave;
[00077] Uw is the wind speed measured at a height of about 19 meters;
[00078] α is 8.0x10-3;
[00079] β is 0.74; and
[00080] g is the acceleration due to gravity.
[00081] The free surface height function f(x) is generated at a point x as follows:

[00082] where for index i > 0,

[00083] and for i < 0,

[00084] In equations (16) and (17), the space wave number for component i is given by
where L is the length of the free surface. The random number N(0,1) can be generated from a Gaussian distribution that has zero mean and a unity variance. As a result, the free surface is formed by adding each wavenumber component that imposes random phase shifts. A Pierson-Moskowitz free surface frozen in time can be computed from equation (16) by using a Fast Fourier Transform for greater computational efficiency.
[00085] Free surface waves are dispersed and in deep water, the frequency and wave number are related by a dispersion ratio given by:

[00086] Equation (18) implies that each spatial harmonic component of the surface can be moved with a defined phase velocity. As a result, in general, longer wavelength free surface waves travel faster than shorter wavelength waves. The combination of Equations (16) and (18) gives a Pierson-Moskowitz free surface that varies over time:

[00087] where t is the instantaneous time. Equation (19) characterizes a one-dimensional rough free surface that moves in the positive x direction and can be used to compute the free surface extent at earlier or later time points.
[00088] In general, a frozen free surface for earlier or later times can be computed as follows. Consider a free surface shape at an instant in time t with wave heights given by f(x,t) the free surface wavenumber F(Ki) spectrum is computed and an arbitrary known dispersion relation
can be used to calculate free surface frozen on previous occasions.
; or later
; of time by:
Free Surface Wave Field Reflectivity Computation
[00089] Returning to Figure 6, in block 606, the free surface wave field reflectivity (ie, the response of a unit point source at the receiver position) can be computed over the extended free surface by using

[00090] where
is the asymptotic form of the first-order Hankel function with n = 0 and 1. The parameters of equation (21) are represented in Figure 13 as follows:
[00091] k is the wave number of the wave field that propagates;
[00092] f(x) is the free surface height described above;
[00093]
is a coordinate position 1302 of a point on the moving scattering surface;
[00094]:
is a vector 1304 from the origin to the moving scattering point 1302;
[00095]
is a vector 1306 from the origin to a fixed receiver position 1308 with spatial coordinates (xr,yr,zr)
[00096]
: is a vector 1310 from the moving spreading point 1302 to the position of the receiver 1408;
[00097] -
is a vector 1312 from the origin to a font position 1314;
[00098]
' is vector 1316 from source position 1314 to receiver 1308;
[00099]
is a vector 1318 from source 1314 to spreading point 1302; and
[000100]
is the slant factor with normal vectors
1320 and
1322 which corresponds to the normal free surface and the direction of the unit vector of the incident field at:
1302 and θ is the angle between the vectors
1320
1322.
[000101] The pressure wave field
' over the entire frequency range can be obtained from pressure measurements on the pressure sensors. However, as explained above with reference to Figure 4C, the recorded particle motion signal can be contaminated with low frequency noise due to vibrations from the trailed cables. As a result, particle motion signals recorded below the threshold frequency may not be useful for determining vertical velocity wave fields at receiver locations along the cable below the threshold frequency.
... Instead, only the vertical velocity wave field
:obtained from measurements of particle motion over a larger frequency range than
normally reliable. A vertical velocity wave field for frequencies below the threshold frequency .-ê.
...., can be calculated from the normal derivative of the pressure wave field
in a process called "low frequency compensation." Wave field decomposition can then be performed to calculate receiver ghost-deleted receiver ghost-eliminated pressure and velocity wave fields. Computing Pressure Wave Field Gradient
[000102] In block 607, the pressure wave field gradient
, is computed. Computing the P pressure wave field gradient
, is determined by the fact that the source depth is greater than or less than a cable depth. Figures 14A to 14B show xz two-dimensional plan views of an exemplary model space associated with a fluid volume. In Figures 14A to 14B, the space includes a representation of a cable 1402 denoted Sr, a source 1404 which produces sound pulses. Cable 1402 and source 1404 may be towed by an exploration seismology vessel (not shown) within the volume of fluid below the free surface that can be imaged as described above to produce an imaged free surface profile represented by curve 1406 and denoted by S0. The imaged free surface S0 includes a free surface profile extension 1408 above the font 1404. Vectors such as vector '
1410 represent spatial source coordinates (xs,ys) vectors such as vector
' 1412 represent coordinates (x,y) in fluid volume and vectors such as vector
1414 represent receiver coordinates. (x,y) along cable 1402. When a depth of source 1404 is less than a depth of cable, as shown in Figure 14A, the expression used to calculate the pressure wave field gradient over the entire frequency range is given by:

[000103] where
is the Fourier transform of the time function from source to source at coordinate location -
and
[000104]
is the free surface wave field reflectivity given by Equation (21).
[000105] To compute free surface wave field reflectivity with a hypothetical source at below the receiver surface and
a knowledge of hypothetical receiver position in
of an extended free surface profile is required. Equation (22) is a first type Fredholm integral equation for the pressure wave field gradient
where the right side of Equation (22) contains only known parameters such as the pressure wave field.
and the free surface wave field reflectivity
Free surface wave field reflectivity with a hypothetical source in
below the receptor surface and a hypothetical receptor in
is computed from the imaged free surface profile. The reflectivity gradient
by the use of numerical gradient techniques. On the other hand, when the source is located at a depth below the Sr cable as shown in Figure 14B, the expression used to calculate the pressure wave field gradient
over the entire frequency range is given by:

[000106] In equation (23), the font function
is not used to calculate
. Note that the solutions of Equations (22) and (23) become unstable when the pressure wave field spectrum has very small values (eg, near receiver ghost discontinuities). These spectral discontinuities generally occur for common towing depths in
As a result, the pressure wave field gradient can be computed for frequencies below the threshold frequency wtH for a time-varying arbitrarily rough free surface. A derivation of equations (22) and (23) is provided in an APPENDIX below.
[000107] Depending on whether the source is at a depth above or below the cable as shown in Figures 14A and 14B, respectively, equations (22) and (23) can be solved numerically for the pressure wave field gradient ' in receiver locations along the cable by using quadrature or expansion methods. For quadrature methods, integrals are approximated by quadrature formulas and the resulting system of algebraic equations is solved. For expansion methods, the solution is approximated by an expansion in terms of basic functions. Compute Normal Velocity Wave Field
[000108] Returning to Figure 6, in block 608, the normal component of the vertical velocity wave field below the threshold frequency
at each receiver location along a cable can be calculated by:

[000109] where p is the density of the fluid; and
[000110]
is a vector normal to a receiver.
[000111] Figure 15 shows a segment of a cable 1502 located below an imaged free surface 1604 in the xz plane of a fluid volume. A normal vector 1508 to cable 1502 in receiver 1506 is given by:

[000112] where is
the pressure wave field gradient in receiver 1606; and
[000113]
is the normal derivative of the pressure wave field P in receiver 1606.
[000114] At this point, the pressure wave field
over the entire frequency range is known from the pressure sensor measurements, the yth vertical velocity wave field
over a frequency range below the low frequency threshold wt is obtained from measurements of particle motion and the vertical velocity wave field below the threshold frequency
is calculated from equation (24). The vertical velocity wave field
over the entire frequency range can be calculated by:

[000115] where
is a weighting function with the property that :
When
when to transition smoothly
from to conform ® increases from to -

[000116]
it's smaller than
to maintain a good signal-to-noise ratio of the hydrophone signal.
[000117] For example, the weighting function
: can be a linear weighting function given by

[000118] In other embodiments, other types of linear weighting functions as well as non-linear weighting functions can be used such as a hanning weighting function. Wave Field Computation with Ghost Removed from Receiver Effect
[000119] In block 609, the recorded pressure wave field and the vertical velocity wave field given by Equation (25) are used to calculate the ghosted pressure wave field from the receiver side (i.e., the upward pressure wave field -• ') over the entire frequency range in the wave-frequency number domain by using:

[000120] The pressure component that goes up in each receiver along the cable is calculated from
per:

[000121] and the vertical velocity wave field with ghost effect eliminated in each receiver in the wavenumber-frequency domain is given by:

[000122] Note that the ghosted pressure wave field eliminated
and the ghosted vertical velocity wave field eliminated
are the pressure and vertical velocity wave fields that go up over the entire frequency range.
[000123] The ghosted pressure wave field can then be transformed from the wave-frequency number domain into the space-frequency domain to obtain:

[000124] followed by a transformation from the space-frequency domain to the space-time domain:

[000125] Similarly, the ghosted vertical velocity wave field with ghost effect can be transformed from the wave-frequency number domain to the space-frequency domain for

[000126] followed by a transformation of the vertical velocity wave field from the space-frequency domain to the space-time domain:

[000127] The Fourier transform inversions of equations (29) to (32) can be performed by using inverse fast Fourier transforms for greater computational efficiency. In short, the receiver with eliminated ghosting or wave fields up
and in receivers they are used to calculate images of underground formations located below the fluid volume.
[000128] Figure 16 shows an illustrative example of a generalized computer system that performs an efficient method for eliminating the source ghosting of a scattered wave field and, therefore, represents a seismic analysis data processing system for the which description is targeted. The internal components of many small, medium, and large computer systems as well as processor-based storage systems can be described in relation to this generalized architecture, although each particular system may have many additional components, subsystems, and the like, parallel systems with architectures similar to this generalized architecture. The computer system contains one or multiple central processing units ("CPUs") 1602 to 1605, one or more electronic memories 1608 interconnected with the CPUs by a CPU/memory subsystem bus 1610 or multiple buses, a first bridge 1612 which interconnects the CPU/Memory Subsystem Bus 1610 with additional buses 1614 and 1616 or other types of high-speed interconnect media that include multiple high-speed serial interconnects. These serial buses or interconnects, in turn, connect the CPUs and memory to specialized processors such as a 1618 graphics processor and with one or more additional 1620 bridges that are interconnected with high-speed serial links or with multiple controllers 1622 to 1627 such as controller 1627 which provides access to several different types of computer readable media such as computer readable media 1628, electronic monitors, input devices and other such components, subcomponents and computing resources. Electronic monitors that include visual display screens, audio speakers and other output interfaces, and input devices that include mice, keyboards, touch screens and other such input interfaces together constitute input and output interfaces that enable the computer system interacting with human users. Computer readable medium 1628 is a data storage device that includes electronic memory, magnetic or optical disk drive, USB drive, flash memory and other such data storage devices. Computer-readable medium 1628 can be used to store machine-readable instructions associated with the computational font ghosting methods described above and can be used to store encoded data during storage operations and by which encoded data can be retrieved during operations. readable by computer systems, data storage systems and peripheral devices. APPENDIX
[000129] A derivation of integral equations that can be used to compute the pressure wave field gradient in a cable receiver
is described with reference to Figures 17A through 17E (See also "Extraction of the normal component of particle velocity from marine pressure data," L. Amundsen et al., Geophysics, Vol. 60, No. 1, pages 212 through 222 (January to February 1995)).
[000130] Figure 17A shows two-dimensional xz plane views of an exemplary model space 1700 associated with a fluid volume. Space 1700 includes a representation of a cable 1702 denoted Sr, and a source 1704 that produces sound pulses. Cable 1702 and source 1704 may be towed by an exploration seismology vessel (not shown) within the fluid volume below an imaged free surface profile represented by curve 1706 and denoted S0. The imaged free surface S0 includes a free surface profile extension 1708 above the font 1704. Vectors such as the vector ; r 1704 represent the spatial coordinates (xs,ys) of the source, vectors such as vector ' 1710 represent coordinates (x,y) in space 1700 and vectors such as vector < 1712 represent receiver coordinates (xr,yr) located at the along cable 1702. A shaded region 1714 in Figure 17A and subsequent Figures 17B to 17E represents a scattering region below cable 1702. The speed of sound in region 1714 is given by:

[000131] where c0 is a speed of sound in a fluid; and
[000132]
is the index of refraction over region 1714.
[000133] The pressure field generated by source 1704 at point 1710 is
and can be transformed from the space-time domain to the space-frequency domain by using a Fourier transform to obtain a pressure wave field.
at point 1710. In practice, the transformation from the space-time domain to the space-frequency domain can be achieved by using a Fast Fourier Transform for greater computational efficiency.
[000134] The constant density acoustic wave equation that characterizes a pressure wave field
caused by a single source that generates sound impulses such as source 1704 located in the spatial coordinates "r is given by:

[000135] where is the Laplacian;
[000136] ® is the angular frequency;
[000137]
is the three-dimensional Dirac delta function that represents the output pulse of sound from the 1704 source; and
[000138] is the Fourier transform of the time function from source to source at coordinate location
[000139] The acoustic wave equation (A-2) characterizes the propagation of the pressure acoustic wave field ' in the fluid volume where the acoustic pressure waves originate from the impulse sound source 1704. Replacement of the equation (A -1) for a speed of sound in the acoustic wave equation (A-2) and redistribution give:

[000140] Now consider the second Green identity of the vector calculus:

[000141] where B and C are both scalar fields derivable twice in a volume V bounded by a closed surface S with a vector pointing outward ^. Figure 17B shows a volume V in space 1700. The volume V rests within a closed surface S which is the surface cable Sr 1702 and a hemispherical cover SR 1716 of radius
. The volume V surrounds the imaged free surface S0, but does not include the scattering region 1714. Figure 17B shows two examples of outward normal vectors 1718 and 1720. An integral equation for the pressure wave field P and the gradient of the VP pressure wave field can be obtained by setting B = P and n 5, in Green's second identity (A-4) to give:

[000142] where " is restricted to points within volume V. Green's function
' characterizes free surface reflections and is a solution of the acoustic wave equation (21) for a unitary impulse source within V and is given by:

[000143] where k0 =w/c0.; and
[000144] r represents a point in 1700 space and r the location of the source of the Green function. JD I i'1 I
[000145] For example, reflectivity
; given by the equation
can be used as the Green's function " in equations (A-5) and (A-6) (ie,
). Replacing equations (A-3) and (A-6) in equation (A-5) gives:

[000146] By letting the radius R of the hemispherical shell SR go to infinity, SR approaches an infinite hemispherical shell as shown in Figure 17C. The integral surface over the surface S in equation (A-7) reduces to an integral equation over the surface of cable Sr. Also, the first term in equation (A-7) goes to zero when the source of the reflectivity represented by Green's function is outside volume V. Because the scattering region at 1714 is located outside volume V, the third term in equation (A-7) also goes to zero, which leaves:

[000147] Equation (A-8) represents a functional relationship between the pressure wave field .J which is obtained from measurements made on pressure sensors along the Sr cable and the pressure wave field gradient . J on the pressure sensors along the Sr cable. This functional relationship between the pressure wave field and the pressure wave field gradient given by equation (A-8) can be used to compute the pressure wave field gradient of pressure r " when only the pressure wave field r " is known.
[000148] Equation (A-8) can also be applied to solve two cases where the 1704 source rests above the level of the 1702 cable or rests below the level of the 1702 cable. When the 1704 source is located within the volume V at a depth between the free surface and the receiver surface as shown in the exemplary representation of Figure 17D according to the deviation property associated with the delta Dirac function, the integral over the volume V on the left side of equation (A-8) reduces to :

[000149] The terms of equation (A-9) can be redistributed to give:

[000150] Equation (A-10) is a Fredholm integral equation of the first type for the pressure wave field gradient .on where the right side of the equation (A-10) contains only known parameters such as the field of pressure wave P-- j and the Green function c-, ; ■ < । described above. On the other hand, when the source is located outside volume V at a depth below the Sr cable, as shown in Figure 17E, equation (A-7) reduces to:

[000151] for the pressure gradient
In equation (A-11), the font function
- j is not used to determine
.
[000152] Although the present invention has been described in terms of particular embodiments, it is not intended that the invention be limited to those embodiments. Modifications within the spirit of the invention will be apparent to those skilled in the art. For example, numerous different computational processing method deployments that perform ghosting in efficient receiver and source can be designed and developed by using a number of different programming languages and computer platforms, and by a number of different deployment parameters that include structure. control, variables, data structures, modular organization and other such parameters. Computational representations of wave fields, operators and other computational goals can be implemented in different ways. Although the efficient wave field extrapolation method discussed above can be performed on a single two-dimensional sampling, the elimination of ghosting in receiver and source can be performed by using multiple two-dimensional sampling or three-dimensional sampling obtained experimentally for different depths can be performed. from greater depths to shallower depths and can be applied in other contexts.
[000153] It is noted that the above description of the disclosed embodiments is provided to enable any person skilled in the art to make or use the present disclosure. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the generic principles defined herein may be applied to other embodiments without departing from the spirit or scope of the disclosure. Therefore, the present disclosure is not intended to be limited to the embodiments shown in this document, but is to be in the broadest scope consistent with the principles and novel features described herein.

权利要求:
Claims (20)
[0001]
1. Method of determining structural information about an underground formation to be performed by a computer system including one or more processors (1602, 1063, 1604, 1605) and one or more data storage devices (1608), the characterized method by the fact that it comprises: computing (608) vertical velocity wave fields associated with marine seismic receivers (906) based on reflectivity of a time-varying roughness-free surface profile (1028) above marine seismic receivers and for frequencies below a threshold frequency using pressure wave fields measured at the receivers; and computing (609) receiver phantom-eliminated wave fields based on the measured pressure wave fields, measured vertical velocity wave fields at the receivers above the threshold frequency, and the computed vertical velocity wave fields; and generating an image of the underground formation using at least in part the receiver ghost-deleted wave fields, the image showing structural information about the underground formation.
[0002]
2. Method according to claim 1, characterized in that computing (608) the vertical velocity wave fields further comprises: decomposing (603) the measured wave fields into wave fields going up and down ; computing (604) a free surface profile image of the free surface as a function of wave fields going up and down; computing (605) a free surface profile extension (1212) of the free surface profile (1208) assuming a frozen free surface and using a free surface spatial roughness spectrum; and compute (606) free surface wave field reflectivity to the free surface profile and to the free surface profile extent (606).
[0003]
3. The method of claim 2, further comprising: computing (607) a pressure wave field gradient at receiver-coordinated locations below the threshold frequency using the wave field reflectivity of free surface; and computing (608) the vertical velocity wave field below the threshold frequency as a function of the pressure wave field gradient at the receiver coordinate locations.
[0004]
4. Method according to claim 2, characterized in that decomposing (603) the measured wave fields into wave fields that go up and down further comprises decomposing the measured pressure wave fields into fields waveforms that go up and down and decompose the measured vertical velocity wave fields into vertical velocity wave fields that go up and down.
[0005]
5. Method according to claim 1, characterized in that computing (609) the receiver ghost-eliminated wave fields further comprises computing receiver ghost-eliminated pressure wave fields as a function of measured vertical velocity wave fields and computed vertical velocity wave fields below the threshold frequency, and compute ghosted vertical velocity wave fields from the receiver as a function of the ghosted eliminated pressure wave fields.
[0006]
6. Method according to claim 1, characterized in that marine seismic receivers (906) form a data acquisition surface.
[0007]
7. Method according to claim 1, characterized in that the threshold frequency is a higher frequency of a frequency range over which measured vertical velocity wave fields have a low signal-to-noise ratio.
[0008]
8. Computer system for determining structural information about an underground formation, the computer system characterized in that it comprises: one or more processors (1602, 1063, 1604, 1605); one or more data storage devices (1628); and a ghosting elimination routine stored in one or more of the one or more data storage devices (1628) which, when executed by the one or more processors (1602, 1063, 1604, 1605), the routine being aimed at: retrieving from the one or more data storage devices (1628), pressure wave field data measured at marine seismic receivers above a frequency range and vertical velocity wave fields measured at receivers above a frequency range limit frequency (602); compute vertical velocity wave fields associated with marine seismic receivers for frequencies below the threshold frequency using the measured pressure wave fields and reflectivity of a time-varying roughness free surface profile above marine cables; computing (609) receiver ghost-eliminated wave fields over the frequency range based on the measured pressure wave fields, the measured vertical velocity wave fields above the threshold frequency, and the computed vertical velocity wave fields; and generating an image of the underground formation using at least in part the receiver ghost-deleted wave fields, the image showing structural information about the underground formation.
[0009]
9. Computer system according to claim 8, characterized in that computing (608) the vertical velocity wave fields for frequencies below the threshold frequency further comprises: decomposing (603) the measured wave fields into fields waves going up and down; computing (604) a free surface profile image (1026) of the free surface (1028) as a function of the wave fields going up and down; computing (605) a free surface profile extension (1212) of the free surface profile assuming a frozen free surface and using a free surface spatial roughness spectrum; and computing (606) the free surface wave field reflectivity to the free surface profile and the free surface profile extent.
[0010]
10. Computer system according to claim 9, characterized in that it further comprises: computing (607) a pressure wave field gradient at receiver-coordinated locations below the threshold frequency using the field reflectivity of free surface wave; and computing (608) the vertical velocity wave field associated with the receivers (906) for frequencies below the threshold frequency as a function of the pressure wave field gradient at the receiver coordinate locations.
[0011]
11. Computer system according to claim 9, characterized in that decomposing (603) the measured wave fields into wave fields that go up and down further comprises decomposing the measured pressure wave fields into pressure wave fields that go up and down, and decompose the measured vertical velocity wave fields into vertical velocity wave fields that go up and down.
[0012]
12. Computer system according to claim 8, characterized in that computing (609) the phantom-eliminated wave fields of the receiver further comprises computing the phantom-eliminated pressure wave fields of the receiver as a function of the measured vertical velocity wave fields and the computed vertical velocity wave fields below the threshold frequency, and compute ghosted vertical velocity wave fields of receiver as a function of the ghosted eliminated pressure wave fields .
[0013]
13. Computer system according to claim 8, characterized in that marine seismic receivers (812) form a data acquisition surface.
[0014]
14. Computer system according to claim 8, characterized in that the threshold frequency is a higher frequency of a frequency range over which measured vertical velocity wave fields have a low signal-to-noise ratio.
[0015]
15. Computer-readable non-transient medium characterized by the fact that it has machine-readable instructions encoded therein to allow one or more processors (1602, 1063, 1604, 1605) of a computer system to perform the operations of: computing fields of vertical velocity waveforms associated with marine seismic receivers for frequencies below a threshold frequency using pressure wave fields measured at the receivers and reflectivity of a time-varying roughness free surface profile above marine cables; compute phantom-eliminated wave fields from receiver based on measured pressure wave fields, measured vertical velocity wave fields at receivers above threshold frequency, and computed vertical velocity wave fields; and generating an image of the underground formation using at least in part the receiver ghost-deleted wave fields, the image showing structural information about the underground formation.
[0016]
16. The medium of claim 15, characterized in that computing the vertical velocity wave fields further comprises: decomposing the measured wave fields into up and down wave fields (603); computing a free surface profile image of the free surface as a function of wave fields going up and down (604); compute a free surface profile extension of the free surface profile assuming a frozen free surface and using a free surface spatial roughness spectrum (605); and compute free surface wave field reflectivity for the free surface profile and the free surface profile extension (606).
[0017]
17. The medium of claim 16, further comprising: computing a pressure wave field gradient at receiver-coordinated locations below the threshold frequency using free surface wave field reflectivity (607); and computing the vertical velocity wave field below the threshold frequency as a function of the pressure wave field gradient at the receiver coordinate locations.
[0018]
18. The medium of claim 15, characterized in that computing the receiver ghost-eliminated wave fields further comprises computing receiver ghost-eliminated pressure wave fields as a function of the receiver-phantom wave fields of measured vertical velocity waveforms and computed vertical velocity wave fields below the threshold frequency, and compute ghosted vertical velocity wave fields from receiver as a function of the ghosted eliminated pressure wave fields.
[0019]
19. Medium according to claim 15, characterized in that marine seismic receivers form a data acquisition surface.
[0020]
20. Medium according to claim 15, characterized in that the threshold frequency is a higher frequency of a frequency range over which measured vertical velocity wave fields have a low signal-to-noise ratio.

类似技术:

---

公开号 | 公开日 | 专利标题

---

BR102013017551B1|2021-08-03|COMPUTER METHOD AND SYSTEM FOR DETERMINING STRUCTURAL INFORMATION ON AN UNDERGROUND FORMATION AND COMPUTER-READABLE NON TRANSIENT MEDIA

---

US9279898B2|2016-03-08|Methods and systems for correction of streamer-depth bias in marine seismic surveys

---

BR102013013084B1|2021-02-17|method for computing images of an underground formation, method for enabling processors and system for measuring wave fields in fluid volume

---

BR102016024737A2|2017-05-02|maritime surveys conducted with multiple source arrangements

---

US10459097B2|2019-10-29|Methods and systems for extrapolating wavefields

---

AU2016222338B2|2021-07-29|Wavefield interpolation and regularization in imaging of multiple reflection energy

---

US10466377B2|2019-11-05|Methods and systems for deghosting marine seismic wavefields using cost-functional minimization

---

US9921325B2|2018-03-20|Wavefield separation based on a matching operator between sensor responses in multi-component streamers

---

US10444386B2|2019-10-15|Methods and systems that determine a velocity wavefield from a measured pressure wavefield

---

US9964656B2|2018-05-08|Methods and systems to remove particle-motion-sensor noise from vertical-velocity data

---

US9829593B2|2017-11-28|Determination of an impulse response at a subsurface image level

---

US9791580B2|2017-10-17|Methods and systems to separate wavefields using pressure wavefield data

---

US20200124755A1|2020-04-23|Correction of source motion effects in seismic data recorded in a marine survey using a moving source

---

BR102014005997B1|2021-09-08|METHOD FOR GENERATING AN IMAGE OF AN UNDERGROUND FORMATION, SYSTEM FOR GENERATING AN IMAGE OF AN UNDERGROUND FORMATION AND COMPUTER-READABLE MEDIA

---

GB2530410A|2016-03-23|Methods and systems to remove particle-motion-sensor noise from vertical-velocity data

同族专利:

---

公开号 | 公开日

---

EP2685288A3|2015-07-08|

---

SG10201510275TA|2016-01-28|

---

MX2013008073A|2014-01-17|

---

US9442209B2|2016-09-13|

---

AU2013206627B2|2017-02-16|

---

US20140016436A1|2014-01-16|

---

EP2685288A2|2014-01-15|

---

MX340795B|2016-06-28|

---

AU2013206627A1|2014-01-30|

---

BR102013017551A2|2015-06-30|

---

EP2685288B1|2021-04-07|

引用文献:

---

公开号 | 申请日 | 公开日 | 申请人 | 专利标题

---

---

US4752916A|1984-08-28|1988-06-21|Dan Loewenthal|Method and system for removing the effect of the source wavelet from seismic data|

---

CA2170561C|1996-02-28|2001-01-30|Raymond Wood|Gas, fire and earthquake detector|

---

GB9906456D0|1999-03-22|1999-05-12|Geco Prakla Uk Ltd|Method and system for reducing effects of sea surface ghost contamination in seismic data|

---

US20020118602A1|2001-02-27|2002-08-29|Sen Mrinal K.|Angle dependent surface multiple attenuation for two-component marine bottom sensor data|

---

GB2379741B|2001-09-18|2003-11-19|Westerngeco Ltd|Method for reducing the effect of Sea-surface ghost reflections|

---

GB2389183B|2002-05-28|2006-07-26|Westerngeco Ltd|Processing seismic data|

---

GB2397378B|2003-01-15|2005-03-02|Westerngeco Ltd|Method for retrieving local near-surface material information|

---

US7123543B2|2003-07-16|2006-10-17|Pgs Americas, Inc.|Method for seismic exploration utilizing motion sensor and pressure sensor data|

---

US7359283B2|2004-03-03|2008-04-15|Pgs Americas, Inc.|System for combining signals of pressure sensors and particle motion sensors in marine seismic streamers|

---

FR2916540B1|2007-05-25|2009-08-28|Cgg Services Sa|SEISMIC EXPLORATION METHOD FOR THE REMOVAL OF GHOSTS DUE TO WATER SURFACE REFLECTIONS, AND METHOD OF PROCESSING SEISMIC DATA FOR THE SUPRESSION OF THESE GHOSTS|

---

US7835224B2|2008-03-31|2010-11-16|Westerngeco L.L.C.|Reconstructing low frequency data recordings using a spread of shallow and deep streamers|

---

US7957906B2|2008-05-07|2011-06-07|Pgs Geophysical As|Method for attenuating low frequency noise in a dual-sensor seismic streamer|

---

US7872942B2|2008-10-14|2011-01-18|Pgs Geophysical As|Method for imaging a sea-surface reflector from towed dual-sensor streamer data|

---

US8902699B2|2010-03-30|2014-12-02|Pgs Geophysical As|Method for separating up and down propagating pressure and vertical velocity fields from pressure and three-axial motion sensors in towed streamers|

---

US8937848B2|2010-07-27|2015-01-20|Cggveritas Services Sa|Methods and systems to eliminate undesirable variations in time-lapse seismic surveys|

---

US9217805B2|2010-10-01|2015-12-22|Westerngeco L.L.C.|Monitoring the quality of particle motion data during a seismic acquisition|

---

US9229123B2|2011-07-25|2016-01-05|Pgs Geophysical As|Method for handling rough sea and irregular recording conditions in multi-sensor towed streamer data|

---

US10459099B2|2011-09-22|2019-10-29|Cgg Services Sas|Device and method to determine shape of streamer|

---

US8634271B2|2012-01-11|2014-01-21|Cggveritas Services Sa|Variable depth streamer SRME|

---

US9435905B2|2012-04-19|2016-09-06|Cgg Services Sa|Premigration deghosting of seismic data with a bootstrap technique|

---

US20140200820A1|2013-01-11|2014-07-17|Westerngeco L.L.C.|Wavefield extrapolation and imaging using single- or multi-component seismic measurements|

---

US20140379266A1|2013-06-25|2014-12-25|Westerngeco L.L.C.|Processing survey data containing ghost data|

---

US9678235B2|2013-07-01|2017-06-13|Pgs Geophysical As|Variable depth multicomponent sensor streamer|WO2013137974A1|2012-03-12|2013-09-19|Exxonmobil Upstream Research Company|Direct arrival signature estimates|

---

US10224900B2|2013-06-07|2019-03-05|Cgg Service Sas|Systems and methods for de-noising seismic data|

---

WO2015023376A1|2013-08-12|2015-02-19|Exxonmobil Upstream Research Company|Low frequency seismic acquisition using a counter rotating eccentric mass vibrator|

---

US9513392B2|2014-02-18|2016-12-06|Pgs Geophysical As|Estimation of direct arrival signals based on predicted direct arrival signals and measurements|

---

US10598807B2|2014-02-18|2020-03-24|Pgs Geophysical As|Correction of sea surface state|

---

US9903966B2|2014-04-14|2018-02-27|Pgs Geophysical As|Seismic data acquisition|

---

GB2527406B|2014-04-17|2020-08-05|Pgs Geophysical As|Methods and systems to separate wavefields using pressure wavefield data|

---

US9791580B2|2014-04-17|2017-10-17|Pgs Geophysical As|Methods and systems to separate wavefields using pressure wavefield data|

---

US9829593B2|2014-08-14|2017-11-28|Pgs Geophysical As|Determination of an impulse response at a subsurface image level|

---

US10444386B2|2014-08-29|2019-10-15|Pgs Geophysical As|Methods and systems that determine a velocity wavefield from a measured pressure wavefield|

---

US9964656B2|2014-08-29|2018-05-08|Pgs Geophysical As|Methods and systems to remove particle-motion-sensor noise from vertical-velocity data|

---

GB2530410A|2014-08-29|2016-03-23|Pgs Geophysical As|Methods and systems to remove particle-motion-sensor noise from vertical-velocity data|

---

US10073183B2|2014-10-20|2018-09-11|Pgs Geophysical As|Methods and systems that attenuate noise in seismic data|

---

CN104375185A|2014-11-04|2015-02-25|中国石油天然气股份有限公司|Method and device for removing surface waves from seismic records|

---

US10345469B2|2016-01-29|2019-07-09|Cgg Services Sas|Device and method for correcting seismic data for variable air-water interface|

---

WO2017155731A1|2016-03-11|2017-09-14|Downunder Geosolutions Pty Ltd.|Method for determining free surface reflectivity for seismic data processing|

---

EP3433642A1|2016-03-24|2019-01-30|Saudi Arabian Oil Company|Simultaneous wavefield reconstruction and receiver deghosting of seismic streamer data using an l1 inversion|

---

US10267936B2|2016-04-19|2019-04-23|Pgs Geophysical As|Estimating an earth response|

---

CN105866838B|2016-05-17|2018-10-16|中国石油天然气股份有限公司|A kind of seismic data low-frequency information compensation method and device|

---

US10871586B2|2017-05-17|2020-12-22|Cgg Services Sas|Device and method for multi-shot wavefield reconstruction|

---

US11105945B2|2018-02-21|2021-08-31|Pgs Geophysical As|Processes and systems that attenuate source signatures and free-surface effects in recorded seismic data|

---

CN109143354B|2018-08-22|2020-03-10|中国石油天然气集团有限公司|Method and device for decomposing seismic waveform characteristics|

---

US10996361B2|2018-09-07|2021-05-04|Saudi Arabian Oil Company|Adaptive receiver deghosting for seismic streamer|

法律状态:
2015-06-30| B03A| Publication of a patent application or of a certificate of addition of invention [chapter 3.1 patent gazette]|

---

2018-12-04| B06F| Objections, documents and/or translations needed after an examination request according [chapter 6.6 patent gazette]|

---

2020-06-09| B06U| Preliminary requirement: requests with searches performed by other patent offices: procedure suspended [chapter 6.21 patent gazette]|

---

2021-06-22| B09A| Decision: intention to grant [chapter 9.1 patent gazette]|

---

2021-08-03| B16A| Patent or certificate of addition of invention granted|Free format text: PRAZO DE VALIDADE: 20 (VINTE) ANOS CONTADOS A PARTIR DE 09/07/2013, OBSERVADAS AS CONDICOES LEGAIS. |

优先权:

---

申请号 | 申请日 | 专利标题

---

US13/545,609|2012-07-10|

---

US13/545,609|US9442209B2|2012-07-10|2012-07-10|Methods and systems for reconstruction of low frequency particle velocity wavefields and deghosting of seismic streamer data|

---